***One sector gravity

clear


********************************
use Gravity_V202211, clear
rename tradeflow_baci value
*******************************************

drop if gdp_d==.
drop if gdp_o==.

drop if gdpcap_o==.
drop if gdpcap_d==.

drop if iso3_o==iso3_d

drop if pop_o<1000
drop if pop_d<1000


sort iso3_o iso3_d year


drop if year<1995
drop if year>2012

sort iso3_o iso3_d year

save Gravity_data_PPML, replace




********************************************************
***************Muendler webpage exercise using Chapter 2 data with intra trade




use "Chapter2Application2.dta", clear
**keep if year == 1986 | year == 1990 | year == 1994 | year == 1998 | year == 2002 | year == 2006
drop if year<1995
drop if year>2011
drop pair_id
rename exporter sorc_iso3
rename importer dest_iso3
replace trade = trade/(10^3)
label var trade "Bilateral trade of goods (current US$ billion)"
rename DIST distw
rename CNTG contig
rename LANG comlang
rename CLNY colony
rename RTA rta
gen byte nafta = (sorc=="CAN" & (dest=="MEX" | dest=="USA")) | (sorc=="MEX" & (dest=="CAN" | dest=="USA")) | (sorc=="USA" & (dest=="CAN" | dest=="MEX"))
replace nafta = 0 if year<1994
replace rta = 0 if nafta==1
label var nafta "Indicator. Equal to 1 if i and j in NAFTA"
label var rta "Indicator that is equal to one when i and j are members of same RTA (except NAFTA)"

*** Naive gravity variables
egen production = sum(trade), by(sorc_iso3 year)
label var production "Source country production of goods (current US$ billion)"
egen marketsize = sum(trade), by(dest_iso3 year)
label var marketsize "Destination country market size of goods (current US$ billion)"

*** Variables in logs
foreach var of varlist trade dist production marketsize {  // zero will be missing
  local lbl : variable label `var'
  local newlbl = "Log " + trim(subinstr("`lbl'","of goods (current US$ billion)","",.))
  gen ln_`var' = log(`var')
  label var ln_`var' "`newlbl'"
  }
	
*** Reference country
* Set iso3 to "ZZZ" for reference oountry so fixed effects of reference country last ones created
* Reference country: United States
gen marketsize_ref_bln = marketsize if dest_iso3=="USA"
replace sorc_iso3 = "ZZZ" if sorc_iso3=="USA"
replace dest_iso3 = "ZZZ" if dest_iso3=="USA"
egen marketsize_ref = mean(marketsize_ref_bln)  // spread to all observations
label var marketsize_ref "Reference country (USA) market size of goods (current US$ billion)"
drop marketsize_ref_bln
gen _imports = trade if sorc_iso3~=dest_iso3
egen imports_dst = sum(_imports), by(dest_iso3 year)
gen selftrade = marketsize - imports
drop _imports imports
label var selftrade "Self trade of goods (current US$ billion)"

*** Source and destination IDs (for PPML-HDFE)
encode sorc_iso3, gen(src_id) label(src_iso3)
encode dest_iso3, gen(dst_id) label(dest_iso3)
qui count if src_id==dst_id & sorc_iso3~=dest_iso3
if `r(N)'>0  disp as err "Source and Destination IDs do not match"


**** GE-PPML Step I: Solve the baseline gravity model ****
*========================================================*

*** Step I.a. Obtain estimates of trade costs and trade elasticities baseline indexes
***           Estimate the gravity model in the "baseline" scenario with the PPML estimator

* Estimate (PPML-HDFE)
* note: ln_distw contig comlang nafta not identified under pair effects
* if no pair effects: ppmlhdfe trade ln_distw contig comlang wto_both rta nafta nafta nafta, absorb(src#year dst#year, savefe) noconst d cluster(src#dst)
ppmlhdfe trade rta nafta, absorb(src#year dst#year src#dst, savefe) noconst d cluster(src#dst)
estimate store benchmark

* Predict source and destination effects: exp_srcfe_bln, exp_dstfe_bln
qui summ __hdfe2__ if dest_iso3=="ZZZ"
local dst_fe_zzz = `r(mean)'
gen double exp_srcfe_bln = exp(__hdfe1__ + `dst_fe_zzz' + _b[_cons])
gen double exp_dstfe_bln = exp(__hdfe2__ - `dst_fe_zzz')

* Predict bilateral trade costs from  panel estimation under pair effects
* if no pair effects: gen double tau_sd_bln = exp(_b[ln_distw]*ln_distw + _b[contig]*contig + _b[comlang]*comlang + _b[nafta]*nafta)
gen double tau_sd_bar = exp(__hdfe3__)
replace tau_sd_bar = exp(0) if __hdfe3__==. & sorc_iso3==dest_iso3  // benchmark trade cost of 0 for self trade
gen double tau_sd_bln = tau_sd_bar * exp(_b[rta]*rta + _b[nafta]*nafta)

* Restrict subsequent analysis to 1994 (start year of NAFTA)
keep if year==2000
save Estimates_gravity_Muendler, replace
* Regress pair fixed effects on select predictors to fill missing pair effects
ppmlhdfe tau_sd_bar ln_distw contig comlang nafta if sorc_iso3~=dest_iso3, absorb(src dst, savefe) noconst d cluster(src#dst)
predict tau_sd_norta, mu
replace tau_sd_norta = exp(0) if sorc_iso3==dest_iso3

* Replace missing trade costs with predictions from standard gravity specification
estimate restore benchmark
replace tau_sd_bar = tau_sd_norta if tau_sd_bar==.
replace tau_sd_bln = tau_sd_bar * exp(_b[rta]*rta + _b[nafta]*nafta) if tau_sd_bln==.
gen ln_tau_sd_bln = log(tau_sd_bln)


**Gen iceberg transport costs and then make them relative to internal ones by dividing for each j by d_ii

gen tau = exp(- ln_tau_sd_bln/4)
sort sorc dest
gen id = tau if sorc==dest
bys sorc: egen  domestic_trade_cost = sum(id)
gen taurel = tau/domestic_trade_cost


drop if sorc==dest

gen lX = log(trade)

gen ltau = log(taurel)


************************************
************************************
************************************
****Figure 2a
************************************
************************************

scatter lX ltau,  ytitle(Bilateral trade in 2000 (log)) ytitle(, size(large)) ylabel(, labsize(medium)) xtitle(Trade costs (log)) xtitle(, size(large)) xlabel(, labsize(medium)) legend(off) graphregion(fcolor(white))
graph save Graph "Tradecosts_trade.gph", replace
graph export  "Tradecosts_trade.eps", as(eps)replace
graph export "Tradecosts_trade.pdf", as(pdf) replace


**Create a measure of iceberg trade costrs by group


gen dum_o=1 if sorc=="ZZZ"
replace dum_o=3 if sorc=="CHN"
replace dum_o=2 if dum_o==.

gen dum_d=1 if dest=="ZZZ"
replace dum_d=3 if dest=="CHN"
replace dum_d=2 if dum_d==.

collapse (mean) taurel [iweight = trade], by(dum_o dum_d)
drop if dum_o==dum_d
reshape wide taurel, i(dum_d) j(dum_o)

edit
save TradeCosts, replace



**************Estimate productivity****************

use  Estimates_gravity_Muendler, clear
keep if year==2000
keep sorc year  exp_srcfe_bln
duplicates drop
rename sorc iso3_o
replace iso3_o = "USA" if iso3_o=="ZZZ"
sort iso3_o year
save tempo, replace


use Gravity_data_PPML, clear
keep if year==2000
keep iso3_o year gdp_o gdpcap_o pop_o gdpcap_ppp_o gdp_ppp_o
duplicates drop
sort iso3_o year
merge iso3_o year using tempo
keep if _m==3
drop _m
drop if exp_==.
rename exp_ FE_o
replace FE_o =log(FE_o)


collapse  pop_o gdp_o gdpcap_o gdpcap_ppp_o gdp_ppp_o FE_o, by(iso3_o)

gen w_o=gdpcap_ppp_o
**gdp_o/(pop_o*10^6)
gen w_o_relUS=w_o/w_o[62]
gen FE_o_relUS=FE_o-FE_o[62]
**gen p_o_relUSLZ=((HTS_d/HTS_d[70])*(1/exp(FE_o_relUS)))^0.25
gen p_o=gdp_o/(gdpcap_ppp_o*pop_o*10^6)
gen p_o_relUS=1
**p_o/p_o[62]

sort iso3_o
gen u_o_relUSA=w_o_relUS /p_o_relUS
gen T_o=exp((FE_o_relUS)/4)*u_o_relUSA
gen T_frechet_o=exp(FE_o_relUS)*u_o_relUSA^4
gen gdpcap_o_rel=gdpcap_ppp_o/gdpcap_ppp_o[62]


keep iso3_o pop_o T_frechet_o gdpcap_o_rel u_o gdp_o T_o

gen lT = log(T_f)/4
gen lgdppc = log(gdpcap_o_rel)

************************************
************************************
************************************
****Figure 2b
************************************
************************************
scatter lT lgdppc,  ytitle(Estimated productivity  rel. to US (log)) ytitle(, size(large)) ylabel(, labsize(medium)) xtitle(GDP pc rel. to US (log)) xtitle(, size(large)) xlabel(, labsize(medium)) legend(off) graphregion(fcolor(white))
graph save Graph "T_gdpcap.gph", replace
graph export  "T_gdpcap.eps", as(eps)replace
graph export "T_gdpcap.pdf", as(pdf) replace



gen dum_o=1 if iso3_o=="USA"
replace dum_o=3 if iso3_o=="CHN"
replace dum_o=2 if dum_o==.
collapse (mean) T_f , by(dum)
sort dum
gen Trel = T_f^0.25
replace Trel=Trel/Trel[3]
save T_onesector, replace
edit
